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Dynamical properties of the elliptical stadium billiard, which is a generalization of the stadium 
billiard and a special case of the recently introduced mushroom billiards, are investigated analytically 
and numerically. In dependence on two shape parameters 8 and 7, this system reveals a rich interplay 
of integrable, mixed and fully chaotic behavior. Poincare sections, the box counting method and the 
stability analysis determine the structure of the parameter space and the borders between regions 
with different behavior. Results confirm the existence of a large fully chaotic region surrounding the 
straight line 5 = 1 — 7 corresponding to the Bunimovich circular stadium billiard. Bifurcations due 
to the hour-glass and multidiamond orbits are described. For the quantal elliptical stadium billiard, 
statistical properties of the level spacing fluctuations are examined and compared with classical 
results. 
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I. INTRODUCTION 

An elliptical stadium is a two-parameter planar do- 
main constructed by adding symmetrically two half- 
ellipses to the opposite sides of a rectangle. In the corre- 
sponding elliptical stadium billiard, the point particle is 
moving with constant velocity within this boundary, ex- 
hibiting specular reflections on the walls. This billiard is 
a generalization of the Bunimovich stadium billiard with 
circular arcs[f|, which is fully ergodic for any length of 
the central rectangle. When the semicircles are replaced 
by half-ellipses, a plethora of new dynamical properties 
emerges, ranging from exact integrability, through mixed 
dynamics with a single chaotic component or with a finely 
fragmented phase space, to broad regions of fully devel- 
oped chaos. Thanks to these properties the elliptical sta- 
dium billiard can be considered a paradigmatic example 
of a Hamiltonian dynamical system, which convincingly 
illustrates both gradual and abrupt transitions induced 
by the parameter variation. 

To appreciate these advantages, the system should be 
looked at as a whole. Several descriptions of the el- 
liptical stadium billiard can be found in the literature 
U> H H IE H However, although detailed and math- 
ematically rigorous, they concentrate on some specific 
billiard properties and consider only restricted parame- 
ter regions. Therefore in the present analysis we take 
into account all possible boundary shapes and investi- 
gate the properties of the elliptical stadium billiard in 
the full parameter space spanned by two variables 5 and 
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7. Some of the preliminary results for special cases with 
the present parametrization were described in Q and 
. The elliptical stadium billiard can also be considered 
as a special case of the recently introduced mushroom 
billiards 

Our interest for the elliptical stadium billiard has been 
enhanced by a number of recent experiments and appli- 
cations. Dynamical properties of classical and quantal 
billiards have important consequences for realistic sys- 
tems in optics as well as in the atomic, mesoscopic and 
solid state physics. Here we quote only a few examples. 
Details of the classical and quantal billiard dynamics are 
the necessary tools for designing properties of semicon- 
ducting microlasers and optical devices used in commu- 
nication technologies 0, [^. Measured conductance 
fluctuations in the semiconductor quantum dots of dif- 
ferent shapes are obtained from the wave functions of 
the corresponding quantal billiards [lj, [l5| . The " atom- 
optics" billiards formed by laser beams allow ultracold 
atoms to move within confined regions of desired shape 
jUJ El ■ In the nonimaging optics, the billiard prop- 
erties determine the shapes of the absolutely focusing 
mirrors 19J. As a consequence of the recent technological 
interest in nanowires and similar extended microstruc- 
tures, different problems with the open billiards have 
also been investigatedpfjl l2lj. including the mushroom 
billiards p^ |. Experiments performed with the stadium- 
shaped microwave resonant cavities [23l [24| could be eas- 
ily extended to stadia with elliptical arcs, predictably 
with interesting results. 

In Section II we define the billiard boundary and de- 
scribe its geometrical properties. In Section III the clas- 
sical dynamics of the elliptical stadium billiard is pre- 
sented. The existence and stability of some selected or- 
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FIG. 1: (a) Shape parameters S and 7 for the elliptical sta- 
dium, (b) Meaning of the angles (j>, 0, a and j3. 



FIG. 2: Shapes of the elliptical stadium in dependence on 5 
and 7. Two of them (5 = 7 = 0.50 and 5 = 0.75,7 = 0.25) 
are Bunimovich stadia. 



bits are investigated by means of the Poincare plots and 
orbit diagrams. Properties of the elliptic islands, their 
evolution and the bifurcations induced by the variation 
of the shape parameters are presented. The extent and 
limits of regions in the parameter space with different dy- 
namical behavior are shown. In Section IV further anal- 
ysis of the Poincare sections is presented. The extent of 
the chaotic region of the phase space is numerically esti- 
mated by means of the box-counting method. In Section 
V we consider the quantum-mechanical version of the el- 
liptical stadium billiard, present selected results for the 
level spacing statistics, and compare them with the clas- 
sical results for chaotic fraction of the phase space. Fi- 
nally, in Section VI we discuss the obtained results and 
propose further investigations. 



II. GEOMETRICAL PROPERTIES OF THE 
ELLIPTICAL STADIUM 

The elliptical stadium is a closed planar domain, whose 
boundary in the x — y plane is defined by means of two 
parameters 6 and 7, satisfying the conditions < 6 < 1 
and < 7 < 00. It is symmetrical with respect to the x- 
and y-axis and is described in our parametrization as 



and for 8 < 1 — 7 the two foci are situated at the points 



F 



(3) 



These expressions enhance the importance of the term 
t = 7 2 — (1 — 6) 2 which is negative for 8 < 1 — 7, posi- 
tive for 8 > 1 — 7, and zero for S = 1 — 7 (Bunimovich 
stadium). For 8 < \x\ < 1, the curvature radius is 



R = 



[(1 - 8) 4 + (\x\ - S) 2 ^ 2 - (1 - S) 2 }} 3 / 2 
7(1 -<5) 4 



(4) 



For < I a; I < S the boundary is flat and the curvature 
radius is R = 00. At the endpoints of the horizontal 
and the vertical axis of the ellipse the curvature radius is 
Ri = 7 2 /(l - 8) and R$ = (1 — 8) 2 /j, respectively. It is 
R = 1 — S for the circular arcs of the Bunimovich stadium. 
Fig. 2 shows some typical shapes of the elliptical stadia 
for 8 between and 1 and for 7 between 0.25 and 1.5. 
The parameter 7 can have any value between zero and 
infinity. However, at 7 = the boundary degenerates 
into a line, and the values of 7 greater than 1.5 do not 
introduce any essentially new features. 
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The meaning of parameters is visible from Fig. 1(a). 
The horizontal diameter is normalized to 2. The vertical 
semiaxis of the ellipse is 7, and the height 27 of the bil- 
liard extends from to 00. The horizontal length of the 
central rectangle is 28, and the horizontal semiaxis of the 
ellipse is 1 — 5. 

The limiting boundary subclasses for some specific pa- 
rameter combinations are as follows. For 8 = the shape 
is elliptical, for 6 = 1 it is rectangular. Especially, for 
6 = 7 = 1 it is a square, and for 5 = and 7 = 1 a circle. 
For (5 = 1 — 7 one obtains the Bunimovich stadiapj with 
different lengths of the central rectangle, which define the 
border separating two distinct boundary classes, one for 
6 < 1 — 7 with elongated semiellipses, and the other with 
6 > 1 — 7 and flattened semiellipses. Another interesting 
case is 5 = 7, where the central rectangle is a square, and 
whose properties were briefly discussed in y|. 

The focal points differ for the two classes. For 6 > 1 — 7 
there are four foci at the points 



III. CLASSICAL DYNAMICS OF THE 
ELLIPTICAL STADIUM BILLIARD 

Dynamics of a classical planar billiard can be exam- 
ined by calculating the points on the billiard boundary 
where impacts and elastic reflections of the particle take 
place. Fig. 1(b) shows the meaning of variables de- 
scribing an impact and appearing in the conditions for 
existence and stability of orbits. Symbols 4> and <f>' , re- 
spectively, denote the angles which the directions of the 
incoming and the outcoming path make with the x-axis. 
The normal to the boundary at the impact point T closes 
the angle 8 = (<f> + <fi')/2 with the x-axis. The angle be- 
tween the normal and the incoming (or outcoming) path 
is (3 = (cf>' — (f)/ 2, and a = {tt/2) — (3 is the angle between 
the tangent to the boundary at the impact point T(x, y) 
and the incoming (or outcoming) path. The slope of the 
normal to the boundary at the impact point is given as 
the negative inverse derivative tan 8 = —X/y'. For the 
iterative numerical computation of the impact points it 
is useful to know the relation between the slopes of the 
incident and the outgoing path, expressed as 
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In our numerical computations, two separate sets of data 
were obtained: coordinates x, y of the impact points on 
the billiard boundary, needed for the graphical presen- 
tation of the orbits and for the computations of the 
orbit stability, and coordinates X and V x , suitable for 
the graphical presentation of the Poincare sections. The 
points P(X, V x ) in the Poincare diagrams are obtained 
by plotting as X the x-coordinate of the point S in which 
the rectilinear path segment crosses the x-axis, and as 
V x = cos <f> the projection of the velocity on the x-axis. 
With an additional assumption concerning the horizontal 
diametral orbit (explained in III.B) the Poincare sections 
are thus completely defined. As some orbit segments do 
not cross the x-axis, the number of points in the Poincare 
diagram can be smaller than the number of the impact 
points on the boundary. In our chosen system of units 
the mass and the velocity of the particle are m — 1 and 
V = 1, respectively, hence both X and V* lie in the 
interval [-1,1]. For billiards with noncircular boundary 
segments |sL l25| such variables are computationally more 
convenient than those containing the arc length variable 
suitable for the billiard boundaries with circular arcs|26|. 
Since X and V x are canonically conjugated variables and 
our billiard is a Hamiltonian system which reduces to 
the collision-to-collision symplectic twist map, the phase 
space and the corresponding Poincare sections are area 
preserving. 

Depending on the choice of the shape parameters and 
on initial conditions, one obtains three types of or- 
bits: periodic orbits, which define fixed points in the 
Poincare diagram; quasiperiodic orbits which define in- 
variant curves surrounding the fixed points, and chaotic 
orbits which fill densely a part of the phase space. As 
described in [2(|, periodic orbits can be stable (elliptic), 
unstable (hyperbolic) and neutral (parabolic). To dis- 
cern these properties for the elliptical stadium billiard, 
we use the criterion stated in |2^] , by which the stability 
of an orbit is assured if the absolute value of the trace of 
the deviation matrix M is smaller than 2, thus if 



-2 < TrA/ < 2. 



(6) 



The deviation matrix of the closed orbit of period n can 
be written as M = M±2M23..M n i, where the 2x2 ma- 
trix Mik for two subsequent impact points Tj and T&, 
connected by a rectilinear path (the chord) of the length 
Pik, is 



M lk = 



Pik 
RiRx 



+ 



Ri sin afc 



sin ai sin Oik 



+ 



sin ctfc 
-Ri 



Sill OCk 

sin ai 



Pik 



Rk sin ai 



(7) 

To examine the properties of the elliptical stadium bil- 
liard in the full parameter space, one may proceed by 
scanning the complete array of 6 and 7 values. However, 
it is often sufficient to make use of a restricted, conve- 
niently defined, one-parameter subspace. Depending on 



FIG. 3: Poincare plots for 8 < 1 — 7, obtained by plotting the 
pairs of variables — 1 < X < 1 and — 1 < V* < 1 for 8 = 0.19 
and various values of 7. 



FIG. 4: Some typical orbits appearing for 8 < 1 — 7. Orbits 
in Figs, (a), (c), (d), (f) and (h) are pantographic orbits with 
n = 4, 3, 2, and 1, respectively. 



the property considered, we choose among the following 
possibilities: constant 6 and changing 7, constant 7 and 
changing S as in or S = 7 as in [jj. 



A. Billiards with 8 < 1 — 7 

This subfamily of elliptical stadia has elongated semiel- 
lipses. The corresponding billiards have been investi- 
gated in US H 0, with the principal aim to estab- 
lish the limiting shapes beyond which the billiard is fully 
chaotic. Parameters a and h used there to describe the 
boundary are connected with our parameters 6 and 7 
through relations 



1 = 



6 = 



h 



(8) 



The results reported in 0, suggest that in the pa- 
rameter space there exists a lower limit above which the 
billiard is chaotic, as consequence of the existence of the 
stable pantographic orbits. In Fig. 3 results of our nu- 
merical calculations of the Poincare sections are shown 
for S = 0.19 and different values of 7 < 1 — 5. This exam- 
ple conveniently illustrates the behavior typical for this 
region of the parameter space. Numerous elliptic islands 
arising from the stable pantographic orbits can be recog- 
nized. The islands due to some other types of orbits are 
also visible. Typical orbits contributing to these pictures 
are shown in Fig. 4, for S = 0.19 and increasing values of 
7. The lowest two pantographic orbits, the bow-tie orbit 
(n = 1) and the candy orbit (n = 2), are shown in Figs. 
4(f) and 4(h), respectively. Three higher period panto- 
graphic orbits are seen in Figures 4(a), 4(c) and 4(d). 
Some other types of orbits are depicted in Figs. 4(b), 
4(e) and 4(g). 

In 0, an earlier conjecture of Donnayj^, stating that 
the lower limit of chaos is set by relations 1 < a < 



FIG. 5: (a) Diagram of the two-dimensional parameter space 
(7, 8). Thick lines denote the border of the fully chaotic 
region including parts C, D and E. (b) Enlarged part of (a) 
showing the emergence of multidiamond orbits of higher n. 
Diamonds, circles, triangles, stars and crosses correspond to 
the cases illustrated on Figs. 3, 6, 8, 9, and 10, respectively. 
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and h > 2a 2 \Ja 2 — 1, was investigated. In 
our parametrization this region is delimited by functions 



FIG. 6: Poincare plots for 8 = 0.19 and various 7, for 7 > 
1-5. 



1 



V4- 2V2' 



7 = 




1 -<j 



-1, 



(9) 



(10) 



7=1 -<5. (11) 

In order to visualize the segmentation of the parameter 
space into regions of different dynamical behavior, in Fig. 
5(a) we plot the pairs of shape parameters as points in the 
7 - 6 diagram. Possible parameter choices are situated 
within an infinitely long horizontal band of height 1. The 
tilted line connecting the points (0,1) and (1,0) holds for 
Bunimovich stadia. The points below and above this line, 
respectively, are the elongated and the flattened elliptical 
stadia. The region delimited by functions (9), (10) and 
(11) is visible as a narow quasi-triangular area denoted 
by D in Fig. 5(a). The obtuse angle of this quasi-triangle 
is situated at the point 7 = 0.48711, 8 = 0.47276. In a re- 
cent paper Del Magno and Markarian pj proved exactly 
that for this region of the parameter space the ellipti- 
cal stadium billiard is fully chaotic (ergodic with K- and 
Bernoulli properties). However, there are strong indica- 
tions that the lowest limit of chaoticity is not consistent 
with this line. It can be probably identified with the limit 
derived in separating regions B and C in Fig. 5(a) 
and resulting from the onset of the stable pantographic 
orbits. Whereas the lower limit of the region D consists 
of two parts, this new limit is made of an infinite number 
of shorter segments, corresponding to all possible panto- 
graphic indices n 6] . Our numerical calculations based on 
the box-counting method (described in the next section) 
confirm the chaoticity of the region C and identify the 
border between regions B and C in the parameter space 
as the lower limit of chaos. 

Another conspicuous feature of billiards with 8 < 1 — 7 
is the stickiness of certain orbits and consequent fragmen- 
tation of the chaotic part of the phase space into two or 
more separate sections. This occurs below a certain lim- 
iting combination of 8 and 7. The exact shape of this 
limit is not obvious. It is probably connected with the 
straight line <5 = 1 — 7V2, pointed out in separating 
in Fig. 5(a) the region A from the region B, which can be 
blamed to the discontinuities in the boundary curvature 
at the points where the flat parts join the elliptical arcs. 
Similar phenomena have been noticed in various billiards 
with circular arcs 1271 l2Sl. 



Billiards with 5 > 1 — 7 



This part of the parameter space, comprising the flat 
half-ellipses, had not been investigated previously, except 



for the brief analysis in Q , later cited in 29j . As an ex- 
ample of the typical behavior, the Poincare sections for 
5 = 0.19 and various 7 > 1 — 8 are shown in Fig. 6. Here 
appear some new types of orbits and we examine their 
existence and stability. These are the diametral orbits 
(horizontal, vertical and tilted) of period two, the hour- 
glass orbit of period four, the diamond orbit of period 
four, as well as the whole family of multidiamond orbits. 
In our further description, we will refer to the impact 
coordinates x and y in the first quadrant (instead of \x\ 
and \y\), with no loss of generality for the obtained re- 
sults, due to the symmetries of the boundary and of the 
considered orbits. 



1. Diametral horizontal and vertical two-bounce orbits 

The horizontal two-bounce orbit obviously exists for 
all combinations of 8 and 7. However, the trace of the 
deviation matrix|26| is equal to 



TrM = 2 



2(^-1 
R 



(12) 



where p — 2 and R = 7 2 /(l — 8), so that the stability 
condition (6) reads 



(13) 



Hence the two-bounce horizontal diametral orbit be- 
comes stable for 



S > 1 - 7 2 



7 



(14) 



Thus, there is a bifurcation at the value S = 1 — 7 2 giving 
birth to the stable diametral orbit. Again we choose 8 = 
0.19 as a typical example, and keeping it constant we 
notice that for this value of 8 the upper limit of chaos 
is 7 = VI - 0.19 = 0.9. In the region 1 - 7 < 8 < 
1 — 7 2 , denoted in Fig. 5(a) by E, there are no periodic 
orbits, and this is the region of full chaoticity. This result 
has been proved in 0], among results for several chaotic 
billiards, and later cited in a discussion of the elliptical 
stadium billiards |2^. (It should be mentioned that for 
this case the definition of the shape parameter in [2jj 
differs from that in Eqs.(8).) 

In the Poincare plot, the invariant points of the hor- 
izontal diametral orbit are defined as (0, ±1), since this 
orbit can be understood as a limiting case of a long and 
thin horizontal hour-glass orbit, which will be discussed 
below in more detail. The Poincare diagrams in Fig. 
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6 show that the corresponding elliptic island arising for 
higher values of 7 above the 7 = yl — S line develops 
into a broad band filled with invariant curves. The cor- 
responding orbit will be shown below in Fig. 11(a). 

Here we mention also the existence of a family of the 
vertical two-bounce non-isolated neutral periodic orbits. 
Their properties are identical to those of the bouncing- 
ball orbits between the flat segments of the boundary in 
the Bunimovich stadium billiard 26] and at present do 
not require our further attention. 



FIG. 7: Multidiamond orbits with n = 1 - 8. 



\/l + 7 2 . The curvature radius at the point x = 1 is 
R = 7 2 /(l — 5). This leads to the trace of the matrix 



TrM = 2 



(20) 



2. Tilted diametral two-bounce orbits 



and to the condition p 2 < Rot 1 + 7 2 < 7 2 /(l — S), thus 
the diamond orbit becomes stable when 



A tilted two-bounce orbit having the impact point 
T(x,y) on the billiard boundary with derivative y' ex- 
ists if the following condition is fulfilled: 



mi 



x = 0. 



This is satisfied if 



(21) 



(15) This limiting curve is shown in Fig. 5(a) as a full line 
dividing the region F from G and the region H from I. 



7 2 <5 



7 2 - (1 -5 f 



(16) 



This point should be on the elliptical part of the bound- 
ary 5 < x < 1, which leads to the condition 



7 



(17) 



However, from the stability condition (13), where the 
chord length is 



p = 2 7 Wl 



6 2 



7 2 - (1 - 5f 



(18) 



arises the restriction 



4- Multidiamond orbits 

The multidiamond orbit of order n is the orbit of period 
2 + 2n, which has two bouncing points at the ends of the 
horizontal axis and In bouncing points on the flat parts 
of the boundary (Fig. 7). Such orbit exists if 



S > 1 - -. 

n 



(22) 



Thanks to the trick known from geometrical optics [Tg| 
which allows the mirroring of the billiard around the flat 
part of the boundary, the calculation of the deviation ma- 
trix for the order n becomes identical to the calculation 
of the matrix for the diamond orbit. Thereby the chord 
length p in (20) should be replaced by 



7 



(19) 



L = np n , 



(23) 



which cannot be fulfilled simultaneously with Eq.(17). 
The conclusion is that no stable tilted two bounce orbit 
can exist. 



3. Diamond orbit 

The diamond orbit of period four shown in Fig. 7(a) 
exists for any parameter choice. It has two bouncing 
points at the ends of the horizontal scmiaxis, and the 
other two on the flat parts on the boundary. To assess 
its stability, one should calculate the deviation matrix 
AI = (Moi-Mio) 2 . The angles contained in the matrix 
are given as sinao = j/p and sinai = 1/p, where p = 



where 



Pn = \ -5+ T 



The trace of the deviation matrix is then 



(24) 



TrM = 2 



2p 2 » 2 
R 



(25) 



with R = 7 2 /(l — S). The resulting condition for the 
stability of the multidiamond orbit is then 
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FIG. 8: Poincare plots for a set of different values S = 7, show- 
ing the appearance of successive stable multidiamond orbits 
of higher order. 



The condition S < x < 1 that this point should lie on the 
elliptical part of the boundary leads to the requirement 



5>1- 



7 > V2(l-S). 



(32) 



5>1- 



1 + 7 2 n 2 ' 



7> 



l-n 2 (l-<5)' 



(26) 



We stress again that the special case n — 1 corresponds 
to the diamond orbit. The limiting curves (26) are plot- 
ted in Fig. 5(a) and are shown enlarged in Fig. 5(b). For 
7 — > 00 the minimal values of 5 after which the multidi- 
amond orbits appear are 



lim 5 = 1 - 



(27) 



The emergence of the stable multidiamond orbits can be 
followed by observing the Poincare sections for a set of 
values <5 = 7 (Fig. 8). The values of 5 for which an orbit 
of new n appears are given as intersections of the straight 
line (5 = 7 with curves defined by (26), and satisfy the 
equation 



n 2 <5 3 - (n z - 1)<5 2 + (5-1 = 0. 



(28) 



Especially, for the diamond orbit (n — 1) this equation 
reads 



(5 3 + (5 - 1 = 



(29) 



and the orbit becomes stable for 5 = 7 = (u — l/u)/ y/3 = 
0.6823278, where 2u 3 = V31 + V27- For the same type 
of boundary the stable two-bounce orbit appeared at S = 
7 = (V5-l)/2 = 0.618034. 



5. The hour-glass orbit 

The hour-glass orbit looks like the bow-tie orbit ro- 
tated by 7r/2. It exists if the coordinates x and y of the 
impact point on the boundary and the derivative y' of 
the boundary at this point satisfy the equation 

2xy' + yy' 2 - y = 0, (30) 
giving as solution the x-coordinate of the point of impact 



x = (5 + 



(l-<5) 2 



7 2 - (1 - sy 



(<J+ V<* 2 + [7 2 -(l-<*) 2 ]) 



(31) 



This limit is shown in Fig. 5(a) as a dashed line sepa- 
rating the region G from I and region F from H. If we 
denote the points with positive y by 1 and the points 
on the negative side by -1, the deviation matrix can be 
calculated as 



M = (MnMi_i) 2 



(33) 



with x given by Eq. (31). The corresponding angle a 
needed in the matrix (7) is given by 



sin a 



7(2; — (5) 



v /(l-<5) 4 + (x-(5) 2 [ 7 2 -(l-(5) 2 ]' 
The chords lengths are 



p = = 2x 



and 



Pi 



(34) 



(35) 



(36) 



where x is given by (31) and y is calculated from (1). The 
curvature radius at this point is obtained by substituting 
(31) into (4). If we define 



P 



(- 



P 



R sin a R sin a R sin a 
the trace of the deviation matrix is 



i?sin( 



TrM = 2 [2(2$+ l) 2 - l] . 



(37) 



(38) 



Whereas the left hand side of the stability condition (6) 
is valid automatically, its right hand side is fulfilled only 
if 



1 < $ < 0. 



(39) 



The left hand side limit of (39) $ = — 1 reduces to the 
existence condition (32) and denotes the case where the 
hour-glass orbit degenerates into the horizontal diametral 
orbit. The right hand side limit $ = is identical to the 
condition 



(40) 



R sin a R sin a R sin a R sin a ' 



7 



FIG. 9: Enlarged parts of the Poincare plots, showing the 
stable island due to the hour-glass orbit, for 8 = 0.19 and 
1.30 < 7 < 1.80. 



FIG. 10: Highly enlarged parts of the Poincare plots, showing 
the stable island due to the hour-glass orbit, for 8 = 0.19 and 
1.87 < 7 < 2.00. 



giving the prescription for the numerical evaluation of the 
limit beyond which the hour-glass orbit becomes unsta- 
ble. This is illustrated on Figs. 9 and 10, where Poincare 
sections are shown for typical examples wit h 6 — 0. 19 
and a set of boundary shapes between 7 = y/2(l — 5) = 
1.2728 (where the stable orbit emerges) and 7 = 1.993 
(where it becomes hyperbolic). The island correspond- 
ing to the hour-glass orbit is shown enlarged, with dif- 
ferent scales in Fig. 9 and Fig. 10. The evolution of the 
hour glass orbit has very interesting properties. At first 
this island appears as a tiny point at the upper end of 
the Poincare diagram, then descends following the value 
= cos0 of its vertical coordinate. When 7 is varied, 
the central island acquires a resonant belt, with period- 
icity which changes from 8 to 7, then to 6, 5 and 4, as 
shown in Fig. 9. Near 7 = 1.79 the orientation of the 
rectangular island changes from tilted to horizontal, af- 
ter which a triangular island appears. Near 7 = 1.90 
(Fig. 10) the triangular island shrinks to infinitely small 
size, and then starts growing again, but with the opposite 
orientation. This is the phenomenon of the "blinking is- 
land" already noticed in [2j| and described in |3(j • With 
further increase of 7 an interesting type of bifurcation 
takes place. The typical islands in the Poincare sections 
corresponding to these changes are visible in Fig. 10, 
and the orbits in Fig. 11. At the value 7 = 1.993, consis- 
tent with the upper stability limit expressed by (40), the 
stable hour-glass orbit in Fig. 11(e) is replaced by an un- 
stable hour-glass orbit, shown in Fig. 11(f) for 7 = 2.02. 
Simultaneously with this, two stable orbits appear. They 
have lost the symmetry of the hour-glass orbit, but in the 
mutual relation to each other retain the left-right mir- 
roring symmetry (Figs. 11(g) and 11(h)). This type of 
behavior, known as the pitchfork bifurcation, has been 
encountered in limagon billiards 31j. 



FIG. 11: (a) Stable orbit into which degenerate both the 
two-bounce horizontal orbit and the hour-glass orbit; (b) The 
hour-glass orbit near the lower limit of stability; (c),(d) Orbits 
due to the resonant belts shown in Fig. 9; (e) Stable hour- 
glass orbit at the upper stability limit; (f) Unstable hour-glass 
orbit beyond the upper stability limit; (g),(h) Two sligthly 
deformed mutually symmetrical orbits, formed in the bifur- 
cation of the hour-glass orbit. 



IV. ANALYSIS OF THE PARAMETER SPACE 
PROPERTIES BY MEANS OF THE BOX 
COUNTING METHOD 

In this section we return to the question of limits within 
which the elliptical stadium billiard is fully chaotic. For 
5 > 1 — 7, our numerical calculations of the billiard dy- 
namics confirmed the result of HEH, stating that the 
billiard is chaotic for 1 — 7 < <5 < 1 — j 2 . For 6 < 1 — 7, 
it was proved in pfl that the billiard is fully chaotic for 
the narrow region determined by (9), (10) and (11) be- 
low the Bunimovich line (5 = 1 — 7. In 0, this band 
was conjectured to be much broader, and its lower limit 
attributed to the emergence of the stable pantographic 
orbits. 

Here we propose to test these limits numerically, by 
means of the box-counting method. We calculate the 
Poincare sections for a chosen pair of shape parameters, 
starting with 5000 randomly chosen sets of initial condi- 
tions and iterating each orbit for 100 intersections with 
the x-axis, thus obtaining 500000 points (A,V^) in the 
Poincare diagram. Having in mind the symmetry of the 
billiard, we plot the absolute values |Vjc|) of the 

computed pairs. Then we divide the obtained diagram 
into a grid of small squares of side 0.01, count the number 
n of squares which have points in them and calculate the 
ratio of this number to the total number N of squares. 
The obtained ratio n/N is denoted by <7 c i ass and deter- 
mines the fraction of the chaotic region in the Poincare 
section. In this way also certain points belonging to in- 
variant curves within the regular islands are included. 
However, their contribution is negligible in comparison 
with the considerably larger chaotic contribution, since 
the average number of points per square is 50 and our 
principal aim is to examine the onset of complete chaos. 
Numerically, the size of the box has been varied until 
the saturation of n/N ratio was achieved, providing a 
consistent and reliable method of analysis of the set ob- 
tained in the numerical experiment. We proceed by cal- 
culating the chaotic fraction for a constant value of 5 
and varying 7. Results are shown in Fig. 12, where the 
classical chaotic fraction q c i ass is plotted against 7 for 
a set of values of S. These pictures show that for each 
S there exist pronounced lower and upper limits of the 
fully chaotic region, characterized by g c iass = 1< When 
one plots these limits in the parameter space diagram in 
Fig. 5(a), one obtains two curves, the upper one within 
the region 5 > 1 — 7, which coincides with Eq. (14) 
and separates the region E from region F, and the lower 
one in the region S < 1 — 7, separating regions B and 
C. Thus our numerical experiments confirmed the exis- 
tence of both limits, the lower one determined in 4] from 
the pantographic orbits, and the upper one suggested in 
[2|. Outside this fully chaotic region (consisting of re- 
gions C, D and E in Fig. 5(a)), the value of the fraction 
<Zciass is characterized by oscillations following the param- 
eter variation, called in [22} the " breathing chaos" . At 
some values there are strong discontinuities, which can 
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FIG. 12: Dependence of qtiass on the shape parameter 7 for 
various values of S, calculated by means of the box-counting 
method. There is a conspicuous difference between the cases 
S — and S — 0.01, typical for stadium billiards. For each 
S 7^ there exists an interval of full chaos, where (/class = 1. 



be traced to bifurcations of various orbits. 



V. THE QUANTUM-MECHANICAL 
ELLIPTICAL STADIUM BILLIARD 

The elliptical stadium billiard can also be considered 
as a two-dimensional quantal system. The Schrodinger 
equation for a particle of mass m moving freely within 
the two-dimensional billiard boundary is identical to the 
Hclmholtz equation for the wave function "F 



h 2 , 

V 2 * = EV, 

2m 



where E is the particle energy. The usual transformation 
to dimensionless variables is equivalent to substituting 



— = 1 and E = k 2 , 

2m 



(42) 



which yields 



dx 2 



dy 2 



*(x,y) = 0. 



(43) 



We solve this equation following the method of Ridell 32] 
and, using the expansion in the spherical Bessel functions 
, obtain the wave functions and energy levels within 
the preselected energy interval. The described method 
yields about 500-3000 levels, depending on the shape. 



Energy level statistics for the elliptical stadium 
billiard 



According to |3J, |35j , statistical properties of the en- 
ergy spectrum reflect the degree of chaoticity in the cor- 
responding classical system. For the elliptical stadium 
billiard it would be interesting to see whether the quan- 
tal calculation confirms the existence of the region of 
fully developed chaos, in the sense of the conjecture of 
Bohigas, Giannoni and Schmit[34j].To assess this corre- 
spondence, we calculate the spectrum and wave func- 
tions for a given pair of shape parameters. The obtained 
spectrum is then unfolded using the method of French 
and Wong[3f3| and the resulting histograms are analyzed 
by means of the Nearest Neighbor Spacing Distribu- 
tion Method (NNSD). To fit the numerically calculated 
distributions, we propose three possibilities: the Brody 



distribution 37] , the Berry-Robnik distribution [38[ , and 
a two-parameter generalization of both Brody and Berry- 
Robnik distributions 



P PR ( s )=e-( 1 -«) s (l-g) 2 



+e-^-^ s q [2(1 -q) + a(u + l) g w+ V] e"" 9 " 



(44) 



which we shall henceforth call the Prosen-Robnik (PR) 
distribution. For u> = 1, P PR (s) reduces to Berry-Robnik 
distribution, and for q = 1 it is identical to Brody dis- 
tribution. It coincides with Wigner distribution if lo = 1 
and q = 1, and with the Poisson distribution whenever 
lu — 0. Here, a is defined as 









a = 








\oj- 



UJ + l 



(45) 



(41) and Q denotes the Incomplete Gamma Function 



Q(a, x) 



T(a) 



(46) 



The derivation of Eq.(44) was based on the factorized 
gap distribution Z PR (s) 



CO - 



,e< + V +1 



(47) 



introduced by Prosen and Robnik[3J|. The gap distri- 
bution is the probability that no level spacing is present 
in the interval between s and s + As. Relation between 
the level spacing distribution P(s) and the gap distribu- 
tion is P(s) = d 2 Z(s)/ds 2 . The function (44) has been 
evaluated in |40j and was used to test the spectra of sev- 
eral types of the billiards of mixed type[25j, |33J. Fitting 
the calculated histograms to this distribution gives two 
parameters qpn and ljpr. According to [3fl| |. the result- 
ing value of <?pr is the variable which corresponds to the 
classical g c iass, the magnitude of the chaotic fraction of 
the phase space. This distribution gives more realistic 
results than the Brody or the Berry-Robnik distribution. 

The Berry-Robnik and Brody distributions have dif- 
ferent behavior for very small spacings: Brody distribu- 
tion vanishes at s = , whereas for the Berry-Robnik 
distribution P BR (0) = 1 — q 2 . Besides, the Berry-Robnik 
parameter q has a well defined physical meaning: quanti- 
tatively it is the fraction of the phase space which is filled 
with chaotic trajectories, whereas the remaining regular 
fraction of the phase space is equal to 1 — q. However, 
the Berry-Robnik distribution is exactly applicable only 
in the semiclassical limit, and we are exploring the com- 
plete spectrum, including the lowest lying levels. 

In [33 Prosen and Robnik argue that distribution (44) 
describes simultaneously transition from semiclassical to 
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FIG. 13: (a)-(k) Histograms of the level density fluctuations 
for spectra of billiards with 5 = 7, fitted to the Prosen-Robnik 
distribution (44); (1) Comparison of the classical chaotic frac- 
tion g c i ass with the values obtained by fitting the histograms 
with Prosen-Robnik, Berry-Robnik and Brody distributions. 

quantal regime and transition from integrability to chaos. 
The two parameters lu and q characterize these two tran- 
sitions, respectively, so that q retains its meaning as the 
chaotic fraction of the phase space, but is applicable also 
to cases far from the semiclassical limit. The diagram 
of (44) in dependence on lu and q has been presented in 

m 

In applying the generalized distribution (44) to the bil- 
liard (1) we hold both parameters lu and q within the 
limits [0,1]. We choose the subset of elliptical stadium 
billiards with 5 = 7, and in Figs. 13(a)-13(k) we show 
the histograms and their fit with P PR . In the region 
where 6 — 7 lies between 0.44 and 0.63 the distribution 
is very close to Wigner and <?pr = 1. According to [34| 
this would correspond to the complete classical chaos. 
In Fig. 13(1) some obtained values of gpR are shown in 
dependence on 5 and compared with the classical result 
obtained with the box-counting method. The values ob- 
tained with the Brody and Berry-Robnik distributions 
are also shown. 

To the question whether the fluctuations of (fciass with 
gradual changes in shape parameter reflect on the quan- 
tal values of q outside the chaotic region, the answer is 
positive, but requires further investigation. Considering 
the importance of bifurcations examined in Section III, 
one can expect the effects similar to those found in the 
oval billiard [|l| . 

VI. DISCUSSION AND CONCLUSIONS 

In this work we have explored classical and quantal dy- 
namics of the elliptical stadium billiards in the full two- 
parameter space, analyzing two distinct groups of these 
billiards, separated from each other by the set of Buni- 
movich billiards. For 5 < 1 — 7 we have confirmed the 
important role of the pantographic orbits in establishing 
the lower bound for chaos. For 6 > 1 — 7 we have ana- 
lyzed in detail the diametral two-bounce orbits and the 



diamond and multidiamond orbits of increasing periodic- 
ity, creating numerous new corresponding elliptic islands 
when the value 6 — 1 is approached. Especially interest- 
ing is the behavior of the hour-glass orbit which is, while 
remaining stable within a large range of increasing values 
of 7, accompanied with a resonant belt with changing 
periodicity. At the value where the orbit becomes un- 
stable, a bifurcation to a pair of stable quasi-hour-glass 
orbits with distorted symmetry appear. The upper limit 
of chaos 5 = 1 — 7 2 is obtained by estimating the sta- 
bility of the horizontal two-bounce orbit, and coincides 
with the limit proposed in . We have thus numerically 
confirmed the existence of a broad region of chaos in the 
parameter space surrounding the straight line belonging 
to the Bunimovich stadia. 

Thereby the most important dynamical phenomena of 
the elliptical stadium billiards in the full parameter space 
are revealed. There remains, however, the open question 
of the stickiness of orbits in the region 5 < 1 — 7, the frag- 
mentation of the phase space and its connection with the 
discontinuities of the boundary curvature, and possible 
importance of the limiting straight line 6 = 1 — 7\/2. One 
possible method to resolve this is to analyze the statistics 
and the phase space properties of the leaking billiards 
|42l |43| . Preliminary investigations show promising re- 
sults, and should be able to explain the foliation and the 
fragmentation of the phase space of the elliptical stadium 
billiard. 

In conclusion, our investigation has revealed dynamical 
properties of a large two-parameter family of stadium-like 
billiards exhibiting a rich variety of integrable, mixed and 
chaotic behavior, in dependence on two shape parameters 
8 and 7, with the special case of 5 = 1 — 7 corresponding 
to the Bunimovich stadium billiard 0. The proposed 
billiard shapes and obtained results could serve as an ad- 
ditional testing ground for the experimental properties of 
semiconducting optical devices and microwave resonant 
cavities. 
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